A pseudo-spectral method for the Kardar-Parisi-Zhang equation. 
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We discuss a numerical scheme to solve the continuum Kardar-Parisi-Zhang equation in generic 
spatial dimensions. It is based on a momentum-space discretization of the continuum equation 
and on a pseudo-spectral approximation of the non-linear term. The method is tested in (1 + 1)- 
and (2 + 1)- dimensions, where it is shown to reproduce the current most reliable estimates of 
the critical exponents based on Restricted Solid-on-Solid simulations. In particular it allows the 
computations of various correlation and structure functions with high degree of numerical accuracy. 
Some deficiencies which are common to all previously used finite-difference schemes are pointed out 
and the usefulness of the present approach in this respect is discussed. 
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I. INTRODUCTION 

Surface growth is a paradigmatic problem of nonequi- 
librium statistical mechanics with widespread potential 
applications such as molecular beam epitaxy, fluid flow 
in porous media or flame fronts 0, |^] . 

For such extended systems, it is a rather intricate ques- 
tion how to connect, in a direct way, microscopic inter- 
actions to the dynamics at mesoscopic or coarse-grained 
scales. Phenomenological models, based on stochastic 
partial differential equations, selected according to sym- 
metry principles and conservation laws, are often capable 
to reproduce various experimental data. The most well- 
known of these models is the one introduced by Kardar, 
Parisi and Zhang (KPZ) ||. The equation introduced by 
these authors engendered an enormous amount of work, 
that addressed the large number of issues related to it 
j|. Yet, many fundamental properties of the KPZ equa- 
tion are still not well understood. Using renormaliza- 
tion group (RG) theory, various authors attempted to 
estimate critical exponents and the upper critical dimen- 
sion [B, @, R p|. The success of this procedure has been 
limited, so far, by the difficulties of RG techniques to 
reach the strong-coupling regime. An alternative route 
to study stochastic partial differential equations, which 
yields an easy and reliable access to critical exponents, 
hinges on the so-called restricted solid-on-solid (RSOS) 
growth models [0, [HI 11 . This approach is based on 
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simple rules for deposition and diffusion of particles on 
a discrete lattice, and it can be implemented in a very 
efficient and fast way on a computer. However, there is 
no a priori guarantee that a RSOS model belongs to the 
universality class of the continuum equation, albeit there 
is a general believe in this sense. Moreover, the corre- 
sponding coupling parameters possess fixed values in the 
RSOS case, a feature that prevents the exploration of the 
entire phase diagram of the KPZ equation. 

Given these premises, numerical integration appears as 
the most direct and definite way to determine the uni- 
versality class of a given continuum equation. So far, 
finite-differences methods have been exploited in the con- 
text of the KPZ equation (lj, |1| 0. In this framework, 
Lam and Shin [ fi"5| have shown the jeopardy of select- 
ing an incorrect discretization in the framework of finite- 
differences algorithms. In this work, we follow a different 
route and propose a numerical scheme based on a pseudo- 
spectral representation. 

Spectral methods, although almost routinely employed 
in fluid mechanics |16|, have not, so far, been used in the 



context of stochastic equations except in Ref. 17 . In 
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this latter work, two of us have introduced in the simple 
(1 + 1) case, a discretization algorithm based on a pseudo- 
spectral scheme which outperforms classical finite differ- 
ences on various respects. The present work generalizes 
the spectral method to the non-trivial (2 + 1) case, where 
no exact results are available and where finite-differences 
methods may be hampered by large discretization effects, 
as we will shortly discuss. Extended numerical simu- 
lations for the (2 + 1) case, are then reported, along 
with comparisons with results based on finite-differences 
schemes. The absence of uncontrolled discretization ef- 
fects also allow us to compute various correlation func- 
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tions and critical exponents in a rather precise and re- 
liable way. Particularly interesting is the computation 
of various height-height correlation functions and of the 
structure factor, which were not previously calculated. 

The plan of the remaining of the paper is as follows. 
In section II, wc address the general features of finite- 
differences discretizations and introduce the new pseudo- 
spectral procedure. Sect. Ill contains definitions of the 



various quantities employed in this work as well as their 
pertinent scaling forms. Results are presented in Sect. IV 
for the (1 + l)-dimensional and in Sect. |v| for the (2 + 1)- 
dimensional cases, respectively. Conclusions and per- 



spectives are included in Sect. VI 



II. THE DISCRETIZED KPZ EQUATION 

We consider a d— dimensional substrate and denote with 
x the d— dimensional vector locating a point on it. The 
surface, which grows on this substrate, is described at 
each time t by its height ft(x,i). In the continuum ap- 
proximation, which disregards overhangs, this dynamics 
generically satisfies a stochastic partial differential equa- 
tion 



d t h = T{h) + r h 



(1) 



where J-{h) is a functional containing various derivatives 
of ft(x, £), and r](x, t) includes all the fluctuations pro- 
duced by interactions among the unresolved microscopic 
degrees of freedom. This noise term clearly influences the 
dynamics at mesoscopic scales and therefore the global 
properties of these coarse-grained surfaces, in particu- 
lar their rough or super-rough nature. In this work, we 
mainly focus on the case of white noise of zero average 
and amplitude D which is delta-correlated in time as well 
as in space 



(»7(x,t)> 







2D5 d (x - x') 5{t - t'). (2) 



In equations (g) and below, the symbol (...) represents 
ensemble averages over different realizations of the noise. 

Constraints based on symmetry principles are helpful 
in reducing the functional T(h) to a few standard pre- 
scribed forms, according to the type of growth process 
it is expected to model. The KPZ equation ||, which is 
among the most common ones, reads 



d t h = vV 2 h- 



n 



(3) 



where v and A denote coupling parameters for the linear 
and non-linear terms respectively. This equation con- 
tains, in addition to a linear diffusion term, a non-linear 



term which takes into account the local growth normal to 
the interface [|, fjj. In the following, fields are assumed 
to be periodic of characteristic length L = V 1 ^ with 
respect to any spatial directions. 

A. Finite difference discretizations 

In finite-differences methods, one discretizes space by 
defining points xj = jAa; (j being a set of d integer in- 
dices) on a cubic lattice of elementary size Ax. In the 
framework of a one step Euler time discretization, the 
KPZ time evolution reads 



ft(xj,t + At) = h(xj,t) + At 

+At 7/(xj,i), 

where one sets 



uV 2 h+^- (Vft) 



r?(x j; t) 



2AD 



(Ax) d At 



ait). 



(4) 



(5) 



Here the factor £ is a noise of zero average and correlation 



(6) 



and it is taken from a uniform distribution between 
and 1/2. The prefactor 



24D 



1/2 

ensures that the noise 



(Ax) d At 

has a second moment identical to that of the Gaussian 
noise integrated over a time interval At Note that 
the use of a uniform distribution makes no difference and 
speeds up the computation as it is well accepted 
the operators containing spatial derivatives in 
must now introduce their finite-differences approxima- 
tions. The Laplacian operator reads 

1 d 

V 2 ft(x) = — J2[h{x + Ax e M ,t) + 




ft(x - Ax e M , t) - 2ft (x, t)}, 



(7) 



Where e M stands for the basis vector along the fi-th di- 
rection. For the non-linear term, different possibilities 
exist which have been argued to lead to different results 
]l5| , p7| . We restrict ourselves to the two following cases: 
the standard one 



/i(x - Ax e^t)], 



(8) 



henceforth referred to as (ST), and the one proposed by 
Lam and Shin jjjsl, henceforth referred to as (LS) jljj 
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(V/i) 2 (x) = — ^{[ft(x + Ax f) - /i(x, i)] 2 + [Mx, *) - " Ai «)f (9) 
+ [h(x + Ax e^t) - h(x,t)] [h(x,t) - h(x - Ax e^t)}}. 



It is worth noting here that the LS discretization does not 
violate a fluctuation-dissipation theorem which is pecu- 
liar of the (1+1)- dimensions, unlike the ST discretization 



B. Pseudo-spectral discretization 

In spectral discretization, the continuum periodic field 
h(x,t) is first expanded in Fourier modes /i q (i) 



h(x,t) = I^T/^e^, 



where 



d d x h(x,t) 



-iqx 



(10) 



(ii) 



and q = {1ixn\l ' L, . . . , 2Ttnd/ L) is a wave vector defined 
by d integers , \i = 1 , . . . , d. Similarly, one expands 
the noise term rj(x,t) and thereafter applies the Fourier 
transform to the continuum equation (|3j). An infinite 
system of coupled Langevin equations is thus obtained 



dt 



-vq 2 hq(t) 



— ]T(k ■ kOMt)M*)<5 k+ k', q + Ut)- (12) 



k,k' 



The Fourier modes rjq satisfy, according to equations (|^), 
the correlation: 



(Va(t)w(t')) = 2VD8^S{t-t' 



(13) 



The spectral discretization then consists in projecting 
the infinite system ( p"2| ) on the vector space of periodic 
functions of period L with only a finite number of non- 
vanishing Fourier modes. To be more specific, it is as- 
sumed that wavevectors q, k and k' in (|l2) belong to a set 
S indexed by integers |n M | < N/2 for [i = 1, . . . , d. In this 
approximation, the dynamical equations conserve their 
original form with finite sums replacing the original infi- 
nite ones. The noise term and the spatial derivatives are 
discretized in a similar way. This approach is thus more 
consistent when compared to finite-differences. Further- 
more, it provides exact results for the Edward Wilkin- 
son equation (EW), which can be read off from the KPZ 
equation (||) with A = 0. In this case the KPZ equa- 
tion is reduced to a set of (N + l) d uncoupled complex 



Langevin equations. With our method the regularization 
is performed in Fourier space and its efficiency in repro- 
ducing the features of the continuum equation should be 
compared to that of finite differences at constant number 
of Langevin equations. Note that the lattice spacing is 
such that a = L/N. 

The temporal discretization for the above complex 
equations is performed by an Euler one time step method 
p5] |. In this case, computation of the linear term as well 
as the noise is straightforward at each time step. On the 
other hand, the computation of the non-linear term 



Xq 



= E [Mq-k)]Mt)Vk(i), (14) 



k,q-keS 



which originates from one-half of the square of gradi- 
ent V/i(x, i), necessitates more care. We use a pseudo- 
spectral algorithm which hinges on the following four 
steps: 

(a) Extension of the momentum space in such a way to 

have M + 1 modes (M > N) per direction instead 
of N + 1 modes. The supplementary Fourier modes 
hq are simply set to zero. 

(b) Computation, in this extended momentum space, of 

the spectrum of gradient V/i(x, t) by sheer multi- 
plication q hq. 

(c) Computation by Fourier transform of V/i(x, t) in real 

space at (M + l) d collocations points xj located on 
a d-dimensional hyper-cubic lattice of elementary 
size L/M. This is followed by the straightforward 
calculation of 1/2 [V/i(x, t)] at these same colloca- 
tions points. 

(d) By a further Fourier transform of these (M + l) d 

values of 1/2 [V/i(x, <)] , one gets the correct sums 
Xq for all wavevectors q G S if M is sufficiently 
large. The values found for external supplementary 
Fourier modes are simply discarded. 

The entire procedure can be efficiently implemented us- 
ing a standard Fast-Fourier package (23) . Without a pre- 
ventive extension of the number of Fourier modes for this 
specific computation, the exact non-linear term \q would 
not be obtained for q £ S: this is the well-known alias- 
ing problem p6| , p4| . The minimum (and hence most 
efficient) choice for M turns out to be 3A^/2 [^6|. The 
reasons for this are detailed in Appendix [A]. 
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III. SOME RELEVANT QUANTITIES. 

Kinetically rough surfaces generated by stochastic differ- 
ential equations such as (g) generally possess scale in- 
variant properties. We recall below some characteristic 
quantities in which this behavior becomes manifest as a 
power law, and define the relevant critical exponents. 

The correlation between fluctuations of height between 
any two points x, x + r located at distance r, is given by 
the height-to-height correlation function 



G n (r,t) 



l/n 

d d x(\h( Xl t)-h(x + r,t)\ n ) I (15) 



where n is, in general, any real positive number. It is 
generally assumed that G2 (r, t) has the scaling form 



G 2 (r,t) = r x M 



where the function M(y) is such that 



1 const y 



CO 





(16) 



(17) 



with (3 = xl z - z is the dynamical exponent, and \ is 
the roughness exponent. Correlations are then of the 
form G<i(t) ~ r x in the limit r <C i 1 ^ 2 , whereas they 
follow the asymptotic limit G^^i) ~ tfi for r t x l z . In 
the absence of multiscaling |^5|, [2(| , the same asymptotic 
scaling behaviors should be clearly envisaged for different 
n's. The law ( |l6| ) can be also connected to the scaling 
of the roughness W(t, L). This latter quantity, which 
measures fluctuations amplitudes, is defined as 



W 2 (t,L) = ((h-h L ) ), 



(18) 



where Kl denotes the average height over a volume 
V = L d . For the Family- Vicsek ansatz p7j |, such a global 
scaling reads 



W(t,L) = L X N 



(19) 



where function N(jj) is such that N(y) — » const when 
y — > 00 and N(y) ~ y 13 when y — > 0. The relation 
between W(t, L) and G2(r,t) is discussed in Appendix 
[b|. The exponents are similar to those of (|l7|): quantity 
P is related to the short time dynamics W(t, L) ~ t° and 
X to the asymptotic saturation width. W(t, L) ~ L x . 
Moreover, note that the roughness may be expressed as 
(see Appendix ^|) 



where S(q, t) denotes the structure factor 

1 



S(q,t) 



V 



h q {t)h- q (t) 



(20) 



(21) 



In momentum space, one infers the scaling behavior to 
be (see Appendix |(]) 



S(q,t) = q-( d +^$(qH), 
where the function <i> is such that 



*(y) 



const for y ^> 1, 
for y <C 1. 



d+2 x 



(22) 



(23) 



Consequently S(q,t) ~ q-( d+2x ) for large t and S(q,t) ~ 

t 2(3+d/z for gmaU t 

In our simulations we computed the skewness param- 
eter s, which is the (scaled) third moment of the height 
fluctuations 



s(t) 



[h-h L y 



3/2 ■ 



(24) 



As observed by Krug et al. [|28| for d — 1, this parameter 
is a good measure of the asymmetry of the distribution 
of height fluctuations. 

Finally, we recall that the KPZ equation carries a close 
analogy with turbulence theory. For instance, the equiv- 
alent of a Reynolds number for the KPZ equation can be 
defined as 12911 



DX d 



(25) 



Moreover, exact heuristic arguments show the existence 
of the following two length scales 



h 



In 



( V 



d+3 \V2 



i/(d+i) 



(26) 



(27) 



For a length I smaller than Z; n , the diffusion linear term is 
predominant as in the case of viscous scales observed in 
turbulent flows. Conversely, the nonlinear term is domi- 
nant for length / 3> Z; n , as for the inertial scales in turbu- 
lent flows. In fluid dynamics, the outer scale is given by 
the largest scale available, which is generally provided by 
the geometry of the specific problem. In the KPZ case 
p9[ , it is given by imposing that for scales l in < I < l out , 
the surface is rough. Note that, since the dimensionlcss 
number e is given by 



lout 

111! 



2(d+l)/(d+3) 



(28) 



there exists an equivalent of the inertial range for large 
equivalent Reynolds numbers e> 1. In this regime, the 
typical fluctuation of scale / is of amplitude hi and varies 
with a typical time U given by 



l/(d+3) 



I 2 

AV < 29 > 
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In the context of the KPZ equation, the inertial range is 
represented by the strong-coupling regime. This corre- 
sponds to e > 1. However, the requirement of the exis- 
tence of a non- negligible inertial region, that is Z ou t 3> hm 
yields t s ^ » 1 with f(d) = (d + 3)/2(d + 1). Since 
f(d) < 1 for d > 1 then e-^ d ) is a slowly increasing func- 
tion of e. The above remark explains the difficulties of 
the numerical simulations to reach unambiguous results, 
due to the long transients which arc present in (2 + 1) 
and higher dimensions, as we shall further discuss below. 



V. CASE (2 + 1)- DIMENSIONS 

Let us consider the most physically relevant case i.e. that 
of bi-dimcnsional substrate d = 2. To this purpose, re- 
sults for various finite differences discretizations and for 
the pseudospcctral discretization are provided for sizes 
up to L — 256 and L = 512 respectively. Computations 
are averaged over a number of configurations typically of 
the order of 50 — 100. The values of v and D are both 
set to 1 throughout the rest of this section. 



IV. A TEST CASE: (1 + 1) DIMENSIONS 



Finite difference 



We compare the performance of the various method in 
the (1 + l)-dimensional case. In this instance, the exact 
value for the steady-state roughness of the continuum 
equation is known to be [B8j 



W(L) 



D 
Up 



(30) 



Such a quantity W(L) has been computed for sizes up to 
L = 1024 and parameters v = 1, A = 3, D = 1, and aver- 
aged over 100 different realization of the noise and many 
different steady-state configurations. We have used (a) 
the finite-differences scheme given by Eq. (gh (ST), (b) 
the finite-differences scheme given by Eq. ( [10| ) (LS) and 
(c) the pseudospectral scheme (PS). For finite-differences 
we have set Ax = 1 corresponding to L = N for the 
pseudo-spectral method. A first comparison is reported 
in Fig. where the quantity ip(L) = ^jYlvjDL W(L) is 
plotted for the three cases. We note that error bars refer 
to fluctuations within the 100 different noise configura- 
tions and that the considered three cases are compared 
within identical statistics. 

Unlike the standard discretization, both the pseudo- 
spectral and the Lam-Shin discretizations yields exact 
values for the amplitude (dashed line in Fig. Q), within 
the error bars. However the pseudo-spectral method 
yields much smaller fluctuations compared to the Lam- 
Shin discretization as it can be inferred by comparing 
the error bars for LS (dashed lines) with those of the PS 
method (solid line) which are barely visible in Fig.|l| being 
of the order of the points size. 

As a further support to this result, we compute the 
time evolution of roughness before saturation, which al- 
lows to obtain the critical exponent (3 whose exact value 
is 1/3. In this case the roughness is averaged over 50 
different configurations in all cases. Fig. ^ depicts the 
result. The pseudo-spectral method follows rather ac- 
curately the exact value of the exponent for at least 3 
decades, whereas the LS method slightly overestimates 
that value. For the pseudo-spectral discretization, we 
have also computed (not shown) the roughness exponent 
X with the two different methods illustrated in the next 
section, always finding an excellent agreement with the 
exact value x = 1/2- 



To the best of our knowledge, no papers have so far 
considered a comparison among various finite differences 
discretizations in dimensions d = 2. As previously men- 
tioned, we envisage two relevant finite-differences approx- 
imations (ST) and (LS) of the KPZ equation. In Fig. |[ a 
comparison of the roughness W(t, L) as a function of time 
t is reported for these two methods and two values of the 
non-linear parameter (A = 3 and A = 4.5). From the fig- 
ure, one can appreciate the presence of the three regimes 
(linear regime, KPZ regime and saturation regime) which 
are characteristic of any numerical solution of the KPZ 
equation in (2 + l)-dimcnsions when starting from a flat 
configuration. Both the above values of A yield a non- 
negligible region of KPZ regime. From Eq. (p5|), one has 
e = A 2 and = A 10 / 6 , which for A = 3 give e = 9 

and e^W = 6.240 . . . respectively. For such a moder- 
ate non-linearity (A = 3) the inertial region is already 
present, and it increases as A increases. However, an in- 
stability already noticed in previous simulations fl2] , |Tij} ] 
is present for larger value of A. The curvature present 
in the A = 4.5 case of Fig. stems from this instability, 
and it is probably associated with an overestimation of 
the nonlinearity. 

Our results for the critical exponents arc in good agree- 
ment with previous results reported in the literature for 
finite-difference schemes. A summary of all these results 
is provided in Table | for completeness |50| . 



B. Pseudo-spectral method 

The roughness ( p0[ ) computed by the pseudo-spectral 
method is depicted in Fig. ||for L = 128, and various val- 
ues of A. It has been computed after averaging over many 
realizations of noise (typically of the order of 50 — 100). 
The value A = 3 numerically corresponds to the optimal 

la M W& in which 



choice used in previous works 1 12j 
the strong coupling regime is well displayed (see Fig. ||). 
The power law for times t < L z has been reassessed (see 
Fig. H|) and exponent (3 is evaluated to be j3 = 0.229±0.05 
for A = 3 and 4.5 (L = 512 is the size used to get this 
result). 

The universal function N{t/L z ) = W(t/L z )/L x is 
plotted for A = 3 in Fig. |^. The exponents z and x can 
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be estimated using a recent method devised in Ref. [ p3[ , 
which allows the calculation of the two exponents simul- 
taneously with good accuracy. This method provides a 
more accurate estimate of the roughness exponent \ with 
respect to the more commonly used procedure of comput- 
ing the saturation value of W at different system sizes. 
We obtain x = 0.37 ± 0.02 and z = 1.67 ± 0.05. 

Two alternative procedures to independently calculate 
the roughness exponent is based on the computation of 
the height- height correlation function G2(r) and of the 
structure factor S(q) at stationarity. Its results are re- 
ported in Figs. fj] and ||. Our best values for the rough- 
ness exponent x are X = 0.38 ±0.02 from the correlation 
function and x = 0.40 ± 0.01 from the structure factor. 
Both these values are, within the error bars, compatible 
with previous numerical results on the continuum KPZ 
equation in real space |lj, [l3|, [l4| |3l], |32l a nd with recent 
extensive simulations on RSOS models|llf| . As explained 
in Ref. , the computation of the structure factor yields 
more reliable results of the exponent x with respect to 
the commonly used procedure of extracting it from the 
value of the saturation roughness at various sizes. 

Fig. [| depicts the time evolution of skewness s(t) for 
different sizes L. For A > 0, the quantity s(t) reaches 
an asymptotic value (0.25 — 0.30) indicating a clear-cut 
asymmetry as opposed to the EW case (A = 0). This 
asymmetry in (2 + 1) is intrinsically different from the 
one observed for short times both on the KS (1 + 1) |3f| 
and for the single-step one-dimensional model (2^] , since 
in that instance the asymptotic value of the skewness 
should be zero in both cases, due to the Gaussian nature 
of the probability distribution function. The existence of 
a well defined peak for s(t) suggests that this quantity 
could be a better indicator of the onset of the KPZ regime 
with respect to the roughness. 

It has been suggested |2^, |2(], |3(| , that possible multi- 
scaling behaviors can be inferred by computing the local 
fluctuations G n (l,t) (where L/N < I < L/2, and the 
extreme case in which I equals the lattice spacing L/N 
corresponds to the average values of powers of the sur- 
face gradient) with respect to time. The reason for this 
can be traced back to the fine balance between differ- 
ent terms present in the Langevin equation, necessary 
to get standard scaling. On that basis, one does not 
expect multiscaling in the KPZ equation case (^6|. A 
direct computation confirms this. First we note that dif- 
ferent moments of the stationary state G n (r) depicted 
in Fig. 10 have identical scaling. Next we consider the 
time evolution of G n (l,t) reported in Fig. [□]. This type 
of calculation has been already performed in Ref. |36j| 
and |^5| in the context of discrete models by considering 
the fluctuations of nearest-neighbors points as a function 
of time. Here calculations are performed for points lo- 
cated few lattice spacings away i.e. I is included within 
a circular region 5a < I < 8a j37j. It is apparent that 
different moments behave in an equal way in agreement 
with scaling arguments based on a Flory theory [|g| . 

As a remark, we note that it would be very interest- 



ing to compare the above result with an analogous cal- 
culation for the Kuramoto-Shivashinsky equation (KS), 
since there are both analytical |3^, |3^] and numerical 
evidences that the (KS) can be mapped, at 
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a coarse grained level, onto a KPZ with a large surface 
tension coupling v. 



VI. CONCLUSIONS 

We have presented a comprehensive study of the pseudo- 
spectral method applied to the numerical solution of the 
KPZ equations. The method can be reckoned as an al- 
ternative and highly reliable procedure to the widely ex- 
ploited finite-difference schemes, and its use in the con- 
text of stochastic partial differential equation is new. At 
the price of a moderate increase in the numerical effort, 
the method offers an improved accuracy in the compu- 
tation of the critical exponents, since it does not suffer 
of the major drawbacks which are common to all the 
finite-difference schemes in real space. The reason for 
this is that the functional form of the continuum equation 
is guaranteed to be maintained, the only approximation 
made in the calculation being the truncation to a finite 
rather than infinite number of modes. We have detailed 
how one can carefully deal with the computation of the 
non-linear mode, by using a back and forth transforma- 
tion between real and momentum space which is more 
efficient than a brute force computation of the non-linear 
term directly in momentum space. We have shown how 
finite-difference schemes lead to non-negligible differences 
in the universal behavior in the temporal region which is 
accessible to numerical simulations, and explained why 
the pseudo-spectral method is the one that both theoret- 
ically and practically most closely approaches the contin- 
uum limit. As a non-trivial application of the method, we 
have shown, using extensive simulations and comparing 
the various different schemes, that the currently most 
accepted results for the critical exponents /) and \ in 
(2 + l)-dimensions, as obtained from RSOS simulations, 
are directly accessible. We have also presented different 
ways of computing the roughness exponent x which yield 
rather precise and consistent results. 

While our results are confined to the KPZ equation in 
(1 + 1) and (2 + 1) dimensions, the method is fully general 
and can be extended to higher dimensions and other non- 
linear continuum equation of the Langevin class. 

We believe that the results presented in this work open 
new perspectives in the computation of the universality 
classes of non-equilibrium problems by avoiding the un- 
controlled use of spatial discretizations from the outset. 
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APPENDIX A: THE ALIASING PROBLEM 

In this appendix we discuss the necessity of extending 
the number of Fourier modes to compute the non-linear 
terms ( |l4|) at each time step, and explain why M = 3N/2 
is the optimal choice for the size of the extended space. 
For simplicity, we will treat the d = 1 case, the extension 
to higher dimension being straightforward. 

Consider the function h(x, t) at time t. By definition, 
it contains TV + 1 non vanishing modes (see Fig. 



h{x, t) 



N/2 

E 

k=-N/2 



(Al) 



with qk = 2irk/L. In dimension d = 1, one has S = 
{q k = 2nk/L, N/2 <k< N/2} and V = L. The com- 
putation of the gradient d x h(x,t) or its square at any 
given point, can be carried out correctly by steps (a)-(c) 
without any extension of the momentum space. Instead, 
step (d) is more delicate and it does require the momen- 
tum space extension. To understand this point, let us 
tackle the following more general question. Consider any 
given periodic function of period L 



u(x, t) 



1 



(A2) 



Assume we want to approximate u(x,t) by a truncated 
function u tr only containing 2P + 1 non vanishing modes 



(x,t) 



1 

I 



E <w ^ 



(A3) 



p=-p 



Since such a function u is completely defined if it is 
known for M = 2P+ 1 points, we choose u tr by imposing 
that it takes the same values as u on M collocation points 
Xj = jL/M j = 0, ■ • • , M — 1. What is the error made 
on the Fourier modes or, equivalently, how different are 
u qp and u q r p for -N/2 <p< N/21 

For collocation points Xj — Lj/M j 
the following relation holds 



0.. 



,M — 1, 



for —P < k < P. Going back to the momentum space, 
the true spectrum is thus somewhat modified: this is 
known as the aliasing problem. Note that large wavenum- 
bers are, in principle, more affected than small wavenum- 
bers. 

In the specific case considered, function u(x, i) equals 
1/2 [d x h(x, t)] and only contains a finite number of non 
vanishing modes: to be precise, the 27V + 1 modes such 
that q k = 2-Kk/L with k = —N, ■ ■ ■ , N (see Fig. |l| and 
|l3| ). Observe that we are interested only in the modes 
with |fc| < N/2. If we use for step (d), 

• a truncated function with M = N + 1 modes 
(P = N/2) to approximate 1/2 [d x h(x, t)f, all the 
Fourier components h qk are modified. 



• a truncated function with M — 2N+1 modes (P = 
N) to approximate 1 /2 [d x h(x, t)] , all the non van- 
ishing Fourier components h qk ( k — —N, • • • , N) 
are correct, even those corresponding to spurious 
modes k = —N, • • • , -N/2 - 1 k = N/2 + 1, • • • , N. 
By definition rM — for r ^ hence the in- 



finite sum in (A5) is reduced to the single term 
r = 0. 



• a truncated function with M — 3N/2 + 1 modes 
(P = 3/V/4) to approximate 1/2 [d x h(x, t)] 2 , the 
Fourier modes of the gradient square q p = 2np/L 
are correctly obtained for k = —N/2, • • • , N/2 and 
spurious modes k = -N, • ■ ■ , -N/2-1; k = N/2 + 
1, ■ • ■ , N are not but this is of no significance. Again 
by definition h qk+rM = for r ^ and \k\ < N/2 
hence the infinite sum in ( |A5| ) is reduced to the 
single term r — for |fc| < N/2. 

As a consequence the choice M = iN/2 + 1 (P = 3/V/4) 
for step (d) is the more economical. 



APPENDIX B: CORRELATION FUNCTION AND 
ROUGHNESS 



E M*)^= E E 



"q P + rM 



e i9 ^.(A4) 



p— — P \r— — oo 



By identification of ( |A3| ) and (A4), one gets 

oo 

E ( A5 ) 



By definition one has: 



G 2 2 (r,t) = \\ v rf <1 x([/ 1 (x + r,<)-/ 1 (x,f)] 2 )(Bl) 
Expanding this expression and integrating over r one gets 



r— — oo 
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~jf d d rG\(r,t) = W 2 (t,L) + ^J^ d d x d d r(h(x,t) h(r,t)) 

+ T72 / ddx d<lv (^( x + r ' *)> - T72 / ^ x rfdr (M x + r ' *) ft ( x ' *)> 



(B2) 



For large V we thus find 

W 2 (t,L) = ± J v d d rGl(r,t). 



(B3) 



It is then easy to show that if Eq. (|lq) and Eq 
hold true, then the above relation implies that Eq. (1£) 
is valid. 



APPENDIX C: SCALING OF THE STATIONARY 
STRUCTURE FACTOR 



Using the Fourier expansion, we get: 



h(x + r,i) - h(x,t) = 



gik-(x+r) _ gik-x 



(CI) 



By substituting in G\ (r, t) one gets after few simple ma- 
nipulations 

Gl(r,t) = 

^ (ftk(«) 7U(t)) 2 [1 - cos (k ■ r)] . (C2) 



Note that, in the above equation, the term k = yields 
a vanishing contrib utio n to the sum. Eq. (C2) can then 
be inserted in Eq. (B3), yielding 



W 2 



= ^E(^(*)^-k(*) 



(C3) 



Therefore, since we define the structure factor from the 
relation 



huAt) h k2 (t)) = V6 kl! - k2 S(k,t), (C4) 



we find 



W 



(t,L) = I^S(k,t). 



(C5) 



k#0 



This relation can be checked to be valid directly from the 
definition of W 2 {t, L). On the other hand from Eq. @ 
we also obtain, using Eq. dcT 



(C6) 



Gl(r,t) = -£S(k,*)(l- 



or inverting 



\J d d rGl(r,t)e~* kr = 

s K0 (x>( k '' f )) - s ( k >*)- 



Therefore for k ^ one gets 

S(k,t) = -\ I d d rG 2 2 {r,t)e- 
2 Jv 



ik-r 



If we then assume that 



Gl(r,t) = r 2 * 5 (-l, 



(C7) 



(C8) 



(C9) 



and using the standard angular integral over the d— di- 
mensional angular variables 



dn d e y ) = (2tt) d/2 _! , 



(Ax) 



one easily gets 



S(M) - fc" (d+2x) $(fc z t) 



,(C10) 



(Cll) 



where 



$ (/c 2 t) = 

\ (2n f 2 jf" dXX^g (^) (A) .(C12) 

It is then immediate to see that 

,,, 2 , const for k z t 3> 1, /7-"ii9 , 1 
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TABLES 



TABLE I: Best estimates of critical exponents f3, \ and z in 
d = 2 as reported in the literature. Note that in the nu- 
merical solution of the continuum equation we have taken 
the value of the non-linear coupling parameters correspond- 
ing more closely to the one used in our work. Here KPZ 
stands for continuum KPZ equation with finite-difference dis- 
cretization, whereas FT stands for field-theoretical methods. 
The error bars in the simulations are typically of the order of 
the last reported digit. 
Model (3 x z Reference 



RSOS 


0.245 


0.393 


1.607 


l3 


RSOS 


0.240 






m 


KPZ 


0.25 


0.39 




ij 


KPZ 


0.240 






ij 


KPZ 


0.240 


0.404 




ni 


Flory 


0.25 


0.4 


1.6 | 




FT 


0.25 


0.4 


1.6 


42] 
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FIGURES 
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FIG. 1: Quantity ip(L) = ^I2v/DLW{L) in (1 + 1)- 
dimensions, computed as a function of the size L. W(L) is 
the steady-state roughness as computed using the standard 
finite difference discretization (ST) (with Aa; = 1), the Lam- 
Shin finite difference discretization (LS) (with Aa; = 1), and 
the pseudo-spectral discretizations (PS) with L — N. All 
these discretizations have an identical number of degrees of 
freedom. The number of configurations used in the average 
is also identical in the three cases. The error bars have been 
distinguished for clarity: grey line for ST, dashed line for LS, 
solid line (barely visible) for PS. 




io 2 io° id 

t 



FIG. 3: Roughness W(t,L) obtained from a finite-difference 
scheme in (2 + l)-dimensions, as a function of time t for 
A = 3, 4.5. Here the lateral size is L — 256, the number 
of configurations is 50, and ST and LS stands for Standard 
and Lam-Shin discretizations respectively. The solid line has 
a slope corresponding to j3 = 0.24. 
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FIG. 4: Roughness W(t,L) as a function of time t for various 
values of A as obtained from the pseudo-spectral methods. 
Here the lateral size is L — 128 and the average is over 100 
different configurations. We have schematically indicated the 
typical three regimes found in the simulations. 



FIG. 2: Roughness W(t,L) in the (1 + l)-dimensional case 
as a function of time t for A = 3. Here the lateral size is 
L = 1024, the number of configurations is 50, and ST, LS 
and PS stands for standard, Lam-Shin, and pseudo-spectral 
discretizations respectively. The solid line has a slope corre- 
sponding to P = 0.33. 
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FIG. 5: W(t, L) in (2 + l)-dimensions, as a function of t when 
A = 4.5 and A = 3. The lateral size is L = 512 and the 
number of configurations is 17. The dashed line has a slope 
corresponding to = 0.23. 
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FIG. 6: Collapse plot of the universal function of N(t/L z ) 
when A = 3 and various sizes L. The dashed line corresponds 
to f3 = 0.23. Parameters are identical to Fig. [l|. The obtained 
values for the exponents are x = 0-37 ± 0.02 and z = 1.67 ± 
0.05. 



FIG. 7: Height-height correlation function G2<V) for A = 3 
and sizes L — 128, L — 512. The dashed line corresponds to 
X = 0.4. 




FIG. 8: Structure factor 5(g) for A = 3 and sizes L 
128,512. The dashed line corresponds to \ = 0.4. 
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FIG. 9: Skewness s(t) for A = 3 and sizes L ranging from 
L = 32 to L = 512. 
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FIG. 12: Illustrating the aliasing problem when M = N + 1. 
The arrow line indicates the modes used in the compu- 
tation, (a) the original field h(x,t); (b) the true func- 
tion 1/2 [d x h(x, t)] 2 ; (c) the truncated function associated to 
1/2 [d x h{x, t)\ 2 when M = N + 1 



FIG. 10: Various moments of the height difference distribu- 
tion at stationarity G n (r) for A = 3 and size L = 512. Note 
that the largest value of r corresponds to L/2. 
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FIG. 11: Various moments of the height-height correlation 
functions G„(l,t) for a fixed I ~ 6L/N and A = 3, L = 512. 
The solid line corresponds to j3 = 0.24. 
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FIG. 13: Solving the aliasing problem with M = 3/2N + 1. 
The arrow line indicates the considered modes, (a) The orig- 
inal field h(x,t); (b) the true function 1/2 \d x h{x, t^] 2 ; (c) 
the truncated function associated to 1/2 [d x h(x, t)] when 
M = 3/2iV + l 
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